/** \file    smooth.h
  * \brief   Image smoothing operations.
  * \date    1998-2001
  * \author  Aloudin (Albert) Akhriev
  * \version $Header: /home/cscom1r4/cvsroot/cs/include/imglib/smooth.h,v 1.1 2004/02/04 14:43:24 cscom1r4 Exp $
  */

#ifndef _smooth_h_
#define _smooth_h_

#include "Arr2Dimg.h"
#include "cinline.h"
#include "math/GaussFilter.h"
#include <float.h>
class QImage;

namespace img {

///simultaneous smoothing and/or differentiation
template <class T>
void smooth(const img::image<T> &img, img::image<T> &s, double smooth_sigma);

template <class T>
void diff(const img::image<T> &img, img::image<T> &dx, img::image<T> &dy, 
		  double smooth_sigma);

template <class T>
void diff2(const img::image<T> &img, img::image<T> &dx2, img::image<T> &dxdy, img::image<T> &dy2,
		   double smooth_sigma);

////////////////////////////////////////////////////////////////////////////
// F U N C T I O N S:
////////////////////////////////////////////////////////////////////////////
//
// void ConvolveImageHorz(const YATL::Arr2D<T> & source, const float * const filter,
//                        int size, bool bSmooth, YATL::Arr2D<T> & dest)
//
// void ConvolveImageVert(const YATL::Arr2D<T> & source, const float * const filter,
//                        int size, bool bSmooth, YATL::Arr2D<T> & dest)
//
// void OrdinarySmoothImage(const YATL::Arr2D<T> & source, const float * fx, 
//                          const float * fy, int size, YATL::Arr2D<T> & image)
//
// void GaussSmoothImage(YATL::Arr2D<T> & image, float sigma)
//
// void GaussDSmoothImage( const YATL::Arr2D<T> & source, YATL::Arr2D<T> & dest, float sigma
// 							bool diffX, bool diffY);
//
// void RepeatedlySmoothImage(T * image, int W, int H, int nRepeat)
//
////////////////////////////////////////////////////////////////////////////

const int  MAX_IMAGE_SIZE = 1024; // max. width or height of an image

//--------------------------------------------------------------------------
// Function convolves image with 1-D filter along X-direction. Source
// image can coincide with destination one. Each 1-D filter should be
// represented by its half, including central sample: f[0],f[1],...f[n-1].
// The function proposes that f[-i] == f[+i] for smoothing filter, and
// f[-i] == -f[i] for differentiating one. While smoothing, if 1-D filter
// tries to address pixel outside of a scan line, then the first or last
// scan's pixel is taken instead of absent one.
//--------------------------------------------------------------------------
template<class T>
void ConvolveImageHorz(
  const YATL::Arr2D<T>& source,     // image to be convolved
  const float * const filter,     // filter
  int                 size,       // size of filter
  bool                bSmooth,    // smooth image, otherwise differentiate
  YATL::Arr2D<T>    & dest)       // out: result
{
  int  W = source.Width();
  int  H = source.Height();
  T ** pDst = NULL;
  T    t, row[MAX_IMAGE_SIZE];
  
  assert(W > size && H > size && W <= MAX_IMAGE_SIZE);
  if (!YATL::HasEqualDim(dest, source))
    dest.Resize(W, H);

  pDst = dest.Data();
  for(int y = 0; y < H; y++)
  {
    const T * pRow = source.BeginRow(y);
    for(int x = 0; x < W; x++)
      row[x] = pRow[x];

    for(x = 0; x < W; x++)
    {
      int L = W-1;
      if (bSmooth)
      {
        t = row[x]*filter[0];
        for(int k = 1; k < size; k++)
        {
          int x1 = x-k;
          int x2 = x+k;
          if (x1 < 0) x1 = 0;
          if (x2 > L) x2 = L;
          t += (row[x2]+row[x1])*filter[k];
        }
      }
      else
      {
        t = T();
        for(int k = 1; k < size; k++)
        {
          int x1 = x-k;
          int x2 = x+k;
          if (x1 < 0) x1 = 0;
          if (x2 > L) x2 = L;
          t += (row[x2]-row[x1])*filter[k];
        }
      }
      pDst[y][x] = t;
    }
  } 
}

//--------------------------------------------------------------------------
// Function convolves image with 1-D filter along Y-direction. Source
// image can coincide with destination one. Each 1-D filter should be
// represented by its half, including central sample: f[0],f[1],...f[n-1].
// The function proposes that f[-i] == f[+i] for smoothing filter, and
// f[-i] == -f[i] for differentiating one. While smoothing, if 1-D filter
// tries to address pixel outside of a scan line, then the first or last
// scan's pixel is taken instead of absent one.
//--------------------------------------------------------------------------
template<class T>
void ConvolveImageVert(
  const YATL::Arr2D<T>    & source,     // image to be convolved
  const float * const filter,     // filter
  int                 size,       // size of filter
  bool                bSmooth,    // smooth image, otherwise differentiate
  YATL::Arr2D<T>          & dest)       // out: result
{
  int  W = source.Width();
  int  H = source.Height();
  T ** pDst = NULL;
  T    t, column[MAX_IMAGE_SIZE];
  
  assert(W > size && H > size && H <= MAX_IMAGE_SIZE);
  if (!YATL::HasEqualDim(dest, source))
    dest.Resize(W, H);

  pDst = dest.Data();
  for(int x = 0; x < W; x++)
  {
    const T * pCol = source.BeginColumn(x);
    for(int y = 0; y < H; y++)
    {
      column[y] = (*pCol);
      pCol += W;
    }

    for(y = 0; y < H; y++)
    {
      int L = H-1;
      if (bSmooth)
      {
        t = column[y]*filter[0];
        for(int k = 1; k < size; k++)
        {
          int y1 = y-k;
          int y2 = y+k;
          if (y1 < 0) y1 = 0;
          if (y2 > L) y2 = L;
          t += (column[y2]+column[y1])*filter[k];
        }
      }
      else
      {
        t = T();
        for(int k = 1; k < size; k++)
        {
          int y1 = y-k;
          int y2 = y+k;
          if (y1 < 0) y1 = 0;
          if (y2 > L) y2 = L;
          t += (column[y2]-column[y1])*filter[k];
        }
      }
      pDst[y][x] = t;
    }
  }
}


//--------------------------------------------------------------------------
// Function smoothes image in ordinary way - convolving by smoothing filter.
// 2-D filter should be separable, i.e. it can be represented by independent
// 1-D filters acting along X and Y axes.
//--------------------------------------------------------------------------
template<class T>
void OrdinarySmoothImage(
  const YATL::Arr2D<T> & source,        // image to be smoothed
  const float    * fx,            // x-filter
  const float    * fy,            // y-filter
  int              size,          // size of filters
  YATL::Arr2D<T>       & image)         // out: result of smoothing
{
  ConvolveImageHorz(source, fx, size, true, image);
  ConvolveImageVert(image, fy, size, true, image);
}


//--------------------------------------------------------------------------
// Function smoothes image using Gaussian filter. While smoothing, if filter
// tries to address pixel outside of a scan line, then the first or last
// scan's pixel is taken instead of absent one.
//--------------------------------------------------------------------------
template<class T>
void GaussSmoothImage(
  YATL::Arr2D<T> & image,     // in: source image,  out: smoothed image
  float      sigma)     // sigma of smoothing
{
  if (image.Size() == 0)
    return;
//  sigma = ALIB_LIMIT(sigma, 1.0f, 50.0f);

  int     W = image.Width();
  int     H = image.Height();
  T **    pImg = image.Data();
  double  scale = 0.0;
  double  revSqSigma = (-1.0)/(2.0*sqr(sigma));
  double  step = 0.01;
  double  halfStep = 0.5*step;
  int     filterSize = (sigma<1) ? 4: (int)(4.0*sigma);
  float * filter = new float [ filterSize ];
  T     * buffer = new T [ __max(W,H) + 2*filterSize ];
  T     * pScan = buffer + filterSize;
  T       t;
  int     k, x, y;
  
  // Compute filter.
  for(k = 0; k < filterSize; k++)
  {
    double x1 = k-0.5;
    double x2 = k+0.5-FLT_EPSILON;
    double sum = 0.0;

    for(double x = x1; x < x2; x += step)
    {
      double t = x + halfStep;
      sum += exp(t*t*revSqSigma);
    }
    filter[k] = (float)sum;
    scale += (k == 0) ? sum : (sum+sum);
  }
  scale = 1.0/scale;
  for(k = 0; k < filterSize; k++)
    filter[k] *= (float)scale;

  // Smooth along horizontal scans.
  for(y = 0; y < H; y++)
  {
    int L = W-1;
    for(x = 0; x < W; x++)
      pScan[x] = pImg[y][x];
    for(x = 1; x < filterSize; x++)
    {
      pScan[-x] = pScan[0];
      pScan[L+x] = pScan[L];
    }
    
    for(x = 0; x < W; x++)
    {
      t = pScan[x] * filter[0];
      for(k = 1; k < filterSize; k++)
      {
        t += (pScan[x+k] + pScan[x-k]) * filter[k];
      }
      pImg[y][x] = t;
    }
  }

  // Smooth along vertical scans.
  for(x = 0; x < W; x++)
  {
    int L = H-1;
    for(y = 0; y < H; y++)
      pScan[y] = pImg[y][x];
    for(y = 1; y < filterSize; y++)
    {
      pScan[-y] = pScan[0];
      pScan[L+y] = pScan[L];
    }

    for(y = 0; y < H; y++)
    {
      t = pScan[y] * filter[0];
      for(k = 1; k < filterSize; k++)
      {
        t += (pScan[y+k] + pScan[y-k]) * filter[k];
      }
      pImg[y][x] = t;
    }
  }
  
  delete [] filter;
  delete [] buffer;
}

/** Smooths image with gauss filter.
  * @param diffX, diffY - if true convolves in that direction not with gaussian but with its derivative
  */
template<class T>
void GaussDSmoothImage( const YATL::Arr2D<T> & source, YATL::Arr2D<T> & dest, float sigma,
	 bool diffX, bool diffY )
{
    GaussFilter gaussFilter;
    gaussFilter.Create( 1, sigma );

	ConvolveImageVert( 
		source, 
		diffY ? gaussFilter.DF(0) : gaussFilter.SF(0),
		gaussFilter.Size(),
		!diffY,
		dest
	);

	ConvolveImageHorz( 
		dest, 
		diffX ? gaussFilter.DF(0) : gaussFilter.SF(0),
		gaussFilter.Size(),
		!diffX,
		dest
	);
}

//--------------------------------------------------------------------------
// Function repeatedly smoothes image with (1, 2, 1) mask along each axis.
// Entire image will be smoothed (including border pixels). Sigma of
// equivalent Guassian roughly equals to: sigma = sqrt((7/12)*nRepeat).
//--------------------------------------------------------------------------
template<class T>
void RepeatedlySmoothImage(
  T * image,       // in: source image,  out: smoothed image
  int W,           // image width
  int H,           // image height
  int nRepeat)     // repeat smoothing 'nRepeat' times
{
  T   tmpBuffer[MAX_IMAGE_SIZE+2];
  T * buf = tmpBuffer+1;
  int i, j, k;

  assert(__max(W,H) <= MAX_IMAGE_SIZE);
  for(; nRepeat > 0; nRepeat--)
  {
    for(j = 0; j < H; j++)   // smooth rows
    {
      T * pRow = image + j*W;

      for(i = 0; i < W; i++)
      {
        buf[i] = pRow[i];
      }
      buf[-1] = buf[0];
      buf[W]  = buf[W-1];
      
      for(i = 0; i < W; i++)
      {
        pRow[i] = buf[i-1] + buf[i] + buf[i] + buf[i+1];
      }
    }

    for(i = 0; i < W; i++)   // smooth columns
    {
      T * pCol = image + i;

      for(j = k = 0; j < H; j++, k += W)
      {
        buf[j] = pCol[k];
      }
      buf[-1] = buf[0];
      buf[H]  = buf[H-1];

      for(j = k = 0; j < H; j++, k += W)
      {
        pCol[k] = (buf[j-1] + buf[j] + buf[j] + buf[j+1]) * (1.0f/16.0f);
      }      
    }
  }
}

template <class T>
void smooth(const img::image<T> &img, img::image<T> &s, double smooth_sigma)
{
	YATL::Arr2D<T> fsrc, fdst;
	img::share_data(fsrc, img);

	img::share_data(fdst, s);
	img::GaussDSmoothImage(fsrc, fdst, smooth_sigma, false, false);
}

///simultaneous smoothing & differentiation
template <class T>
void diff(const img::image<T> &img, img::image<T> &dx, img::image<T> &dy, double smooth_sigma)
{
	YATL::Arr2D<T> fsrc, fdst;
	img::share_data(fsrc, img);

	img::share_data( fdst, dx );
	img::GaussDSmoothImage( fsrc, fdst, smooth_sigma, true, false );

	img::share_data( fdst, dy );
	img::GaussDSmoothImage( fsrc, fdst, smooth_sigma, false, true );
}

template <class T>
void diff2(const img::image<T> &img, img::image<T> &dx2, img::image<T> &dxdy, img::image<T> &dy2,
		   double smooth_sigma)
{
    GaussFilter gaussFilter;
    gaussFilter.Create(1, smooth_sigma);

	YATL::Arr2D<T> fsrc, fdx2, fdxdy, fdy2;
	img::share_data(fsrc, img);

	img::share_data(fdx2, dx2);
	ConvolveImageHorz(fsrc,  gaussFilter.DDF(0), gaussFilter.Size(), true, fdx2);
	ConvolveImageVert(fdx2,  gaussFilter.SF(0),  gaussFilter.Size(), true, fdx2);

	img::share_data(fdxdy, dxdy);
	ConvolveImageHorz(fsrc,  gaussFilter.DF(0),  gaussFilter.Size(), false, fdxdy);
	ConvolveImageVert(fdxdy, gaussFilter.DF(0),  gaussFilter.Size(), false, fdxdy);

	img::share_data(fdy2, dy2);
	ConvolveImageHorz(fsrc,  gaussFilter.SF(0),  gaussFilter.Size(), true, fdy2);
	ConvolveImageVert(fdy2,  gaussFilter.DDF(0), gaussFilter.Size(), true, fdy2);
}

} // namespace img

#endif //_smooth_h_
